Theory of equilibrium flux lattices in unconventional superconductors 
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We investigate equilibrium fiux lattice structures in superconductors with unconventional order 
parameters, such as high- Tc cuprates, using a generalized London model with non-local electro- 
dynamics derived from a simple microscopic model. We find a rich phase diagram containing 
triangular, centered rectangular and square lattices with various orientations relative to the ionic 
lattice, as a function of magnetic field and temperature. 



Existence of a mixed state, characterized by a regu- 
lar array of magnetic flux lines penetrating the mate- 
rial, is perhaps one of the most striking properties of 
type II superconductors. The original pioneering work 
of Abrikosov [||, based on the solution of Ginzburg- 
Landau (GL) equations near the upper critical field Hc2^ 
predicted a triangular flux lattice (FL). This prediction 
was subsequently verified by low field magnetic decora- 
tion experiments on a variety of conventional supercon- 
ductors. In some compounds neutron scattering exper- 
iments revealed deviations from perfect triangular lat- 
tices in stronger fields |Q which where attributed to 
anisotropics in the electronic band structure and other 
effects and were modeled by GL theories containing ad- 
ditional higher order derivative terms reflecting the ma- 
terial anisotropics 

One would expect even richer behavior of flux lattices 
in the new class of heavy fermion and copper-oxide su- 
perconductors as these exhibit highly anisotropic elec- 
tronic structures and, very likely, order parameters with 
unconventional symmetries involving nodes in the gap. 
In high-Tc cuprates much of the experimental and the- 
oretical effort has been focused on the sizable region of 
the phase diagram just below Tc{H) in which the vor- 
tex lattice properties are dominated by thermal fluctua- 
tions Q. While understanding the physics of this fluc- 
tuation dominated regime poses an intriguing and diffi- 
cult statistical mechanics problem, investigation of the 
equilibrium vortex lattice structures at low temperatures 
may provide clues about the microscopic mechanism in 
these materials. So far such investigations have been lim- 
ited to YBa2Cu307_5 compound j|J^ , reveahng vortex 
lattices with centered rectangular symmetry and vari- 
ous orientations relative to the ionic lattice. These have 
been modeled by phenomenological GL theories appro- 
priate for anisotropic superconductors, containing addi- 
tional quartic derivative terms ||^ or a mixed gradient 
coupling to an order parameter with different symmetry 
These works found structures in qualitative agree- 
ment with experiment, but their inherent shortcoming is 
the large number of unknown phenomenological param- 
eters and the subsequent lack of predictive power. Also, 
such GL theory has only been solved for vortex lattice 



near Hc2 , which is experimentally inaccessible in cuprates 
away from Tc- [See however Ref. j^] which holds promise 
for full solution at any field.] We have recently for- 
mulated a generalized London model [ p^ which is valid 
in experimentally accessible region of intermediate fields 
Hci < H Hc2- This model is also phenomenological 
and it contains one unknown parameter which controls 
the strength of the symmetry breaking term. With in- 
creasing magnetic field this model predicts a transition 
from triangular to square FL. While no direct experi- 
mental evidence exists in cuprates at present to confirm 
such a prediction, a similar transition has been recently 
observed in a boro-carbide material ErNi2B2C jO| and 
has been described by a similar London model [[12| . 

In this letter we present a microscopic derivation of 
the generalized non-local London model for an uncon- 
ventional superconductor. Based on this model we for- 
mulate, for the first time, quantitative and largely pa- 
rameter free predictions for the behavior of the vortex 
lattice structure as a function of temperature and mag- 
netic field. Our theory is valid in a large part of the H-T 
phase diagram, only restricted by the inherent domain 
of validity of the London model, H ^ Hc2, and T low 
enough that the thermal fluctuations are unimportant. 
The central result of this work is a prediction that the 
FL geometry in unconventional superconductors will dis- 
play a rich and distinctive behavior as a function of field 
and temperature, undergoing a series of transitions and 
crossovers, and eventually attaining a universal limit at 
T = 0. The London free energy in this limit is non- 
analytic and its long wavelength part is fully determined 
by the nodal structure of the gap function. Such behav- 
ior is caused by the low-lying quasiparticle excitations 
within the nodes and thus could never occur in conven- 
tional superconductors with anisotropic band structures. 

In general the relation between the supercurrent j 
and the vector potential A of the magnetic induction 
B = V X A is non-local in a superconductor, reflecting 
the finite spatial extent of the Cooper pair ~ [ P^ - 
Non-local corrections to physical quantities, such as the 
effective penetration depth, will be of order k~^, where 
K = Ao /^o is the GL ratio and Aq is the London penetra- 
tion depth. For strongly type II materials {k ^ 1) such 
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corrections are negligible. Since cuprate superconductors 
fall well within this class (k is in excess of 50 for most) 
local electrodynamics is always used. However, a closer 
examination suggests that this might not be justified in 
all situations, if, as it is widely believed, these materi- 
als exhibit nodes in the gap. In such a case in place of 
the usual coherence length = wf/ttAq one is forced 
to define an angle dependent quantity, ^o(p) = vf/t^^p, 
which diverges along the nodes. Clearly, in the vicinity 
of nodes the condition Ao/Co(p) ^ 1 is no longer satis- 
fied and, in fact, the extreme non-local limit is achieved. 
Non-local corrections therefore cannot be dismissed in 
unconventional superconductors, especially at low tem- 
peratures when quasiparticles selectively populate nodal 
regions. From the above argument it is also clear that 
such corrections will be highly anisotropic and will in 
general break the rotational symmetry of the flow field 
around the vortex, contributing an anisotropic compo- 
nent to the inter-vortex interaction in the mixed state. 
Thus, on very general grounds, one may expect non- 
triangular FL structures in unconventional superconduc- 
tors. We now illustrate this idea by computing the FL 
geometry for the simplest case of a dx2-y2 superconduc- 
tor with cylindrical fermi surface in the clean limit. 

The non-local relation between j and A is conveniently 
written in Fourier space IT^ ] 



jk = -(c/47r)Q(k)Ak. 



(1) 



Here Q(k) is the electromagnetic response tensor which 
can be computed, within the weak coupling theory, by 
generalizing the standard linear response treatment of 
Gorkov equations flj] to an anisotropic gap. We find 




(2) 



where 7k = v^? • k/2, London penetration depth Aq"^ — 
'l'!Te'^VpN{0)/c'^, Matsubara frequencies LUn = 7rT(2n — 1) 
and the angular brackets mean the Fermi surface aver- 
aging. Eq.(y) is valid for arbitrary Fermi surface and 
gap function. For isotropic gap one recovers an expres- 
sion recently derived by Kogan et al. from the Eilen- 
berger theory [ p^ . One may simplify solving the London 
equation by writing it in terms of magnetic induction 
only. Eliminating j from Eq.(|l|) using the Ampere's law 
j = (c/47r)V X B, one obtains 



Bk-kx [Q-i(k)(kxBk)] -0. 



(3) 



For many purposes it is also convenient to write down 
the corresponding London free energy Fj^, such that 
SFl /(5Bk = gives the above London equation: 

Fl = 5][Bk' + (k X Bk)Q-^(k)(k X Bk)]/87r. (4) 



Here A = A(T) is the temperature dependent penetration 
depth (given below) for which it holds that A(0) = Aq. 

One may study FL structure using this formalism pro- 
vided the cores occupy only a small fraction of the total 
volume, i.e. when H <C i?c2 and n ^ 1. To account 
for the topological winding of the phase around the core 
]lO| , p^ it is then necessary to insert source terms pk on 
the right-hand side of Eq. (H) . A commonly used from is 
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where the prefactor insures proper flux quantization (0o 
is the flux quantum and Q is the area of the FL unit 
cell). For a given kernel Q(k) the FL symmetry is then 
determined by minimizing the Gibbs free energy Gl = 
Fl — HB/4tt (where B = (po/n is the average induction). 

At long wavelengths Q(k) can be evaluated by expand- 
ing expression (^) in powers of 7^. The zeroth order term 
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(6) 



is just the temperature dependent penetration depth 
which, at low temperatures, has the well known T-linear 
behavior 0: A'^ « - (21n2)r/Ad] for a d^2_y2 

superconductor with = Ac;(p^ ~Py) ^^id a cylindrical 
Fermi surface. From now on we shall focus on this simple 
case. The leading non-local term is quadratic in k: 
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The expression Qij = 6ijX^^ + Q\j is easily inverted to 

leading order in k: Q^-^ w X'^[Sij — A^Q^-^^J. Substituting 
this into Eq.(Q) and specializing to fields along the z- 
direction we have 

Fl=J2 ^k[l + A'/c' + X^ficik^ + C2klkl)]/8n. (8) 



Here ^ — vp/i^Ad and A^ is assumed to be a tempera- 
ture dependent solution to the appropriate gap equation. 
Dimensionless coefficients ci and C2 are given by 



(a;2+A2)5/2' 
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Fermi surface has been explicitly parameterized by the 
angle 9 between p and the x axis: = (cos 6*, sin 6*) and 
Ap = Ad cos 29. Coefficients ci and C2 depend on tem- 
perature through a dimensionless parameter t = T / Ad- 
From Eq.(|) one can deduce their leading behavior in the 
two limiting cases: for i <C 1 we find 



It is easy to see that in the local limit Qij{k 
SijX^^ one recovers the usual London free energy 
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FIG. 1. Two high symmetry orientations of centered rect- 
angular unit ceU. 



and for i ^ 1 (i.e., near T^. 
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where a = C(5)(l - 2-5)/87r2 = 0.01272. In the above A 
also depends on i, but this will be unimportant for the 
following qualitative discussion. 

The free energy (||) formally coincides with the one de- 
duced previously from phenomenological considerations 
| pO| . An interesting new feature emerging from the 
present microscopic model is a sign reversal of C2 at inter- 
mediate t implied by Eqs. ( |lT| ) and (|o|). At high temper- 
atures C2 is found to be positive, in agreement with pO[ . 
It has been shown that such a term leads to centered rect- 
angular FL structure with principal axes oriented along 
X or y axes of the ionic lattice (cf. Fig. |l]a). The mag- 
nitude of the distortion from a perfect triangular lattice 
(/3 = 60°) is controlled by the magnitude of C2 and grows 
with increasing magnetic field. Eq.(pTl) shows that at 
fixed field this distortion will initially grow with decreas- 
ing temperature. At low temperatures Eq.(]lO|) predicts 
C2 < 0. This will lead to the same centered rectangular 
FL rotated by 45° (cf. Fig. |^). Numerical evaluation of 
Eq.(^ shows that C2 passes through zero at t* ~ 0.19. 
At this temperature the free energy (^) is isotropic and 
the FL will be triangular at all fields. The sign reversal 
of C2 reflects the competition between the two terms of 
different symmetry in W2 and is a unique consequence of 
the gap function having nodes. 

Another consequence of nodes is the fact that, as it can 
be seen from Eg. (p^ , both ci and C2 diverge as 1/t for 
t —^ 0. This divergence signals that the response tensor 
Q(k) is a non- analytic function of k at T = and the 
expansion in powers of 7^ breaks down. Formally this is 
caused by the fact that at T = at the nodal point the 
expression (|^) for Qij (k) contains a term proportional to 
1/7^. At T = the frequency sum in (2|) becomes an 
integral which can be evaluated exactly with the result: 
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where y — 7k /Ap. For small k the dominant con- 
tribution to the angular average comes from the close 
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FIG. 2. Equilibrium angle /3 as a function of reduced tem- 
perature t — T/Ad for various fields. Open symbols mark 
lattice with orientation along x or y direction while solid sym- 
bols mark the lattice rotated by 45° . We use Aq = 1400A and 
K, — 68. Inset: /9 as a function of field at fixed T. 



vicinity of nodes and can be evaluated by linearizing 
Ap around the nodes. One finds that the leading non- 
local contribution is linear in k rather than quadratic. 



For Q, 



Sij Aq ^ 

(1) 



have 



-M(fc>Co) and Qx 



+ Qij', we 

= Qyx = -^(fc<^o)sgn(fcxfciy), where 
A;> = maxdfca;!, \ky\) and fc< = min(|A:a;|, \ ky\). Prefactor 

/i — tt'^ /8^/2 — 0.8723 is exact in the sense that all cor- 
rections to Qij are 0(fc2). The resulting free energy at 
T = is 



F, 



[1 + Ao 



2l2 



M(A2eo)fc>(fc^-A;2)]/8^. (13) 



The non-local term is clearly non-analytic in k. Its func- 
tional form is universal in the sense that it is independent 
of the Fermi surface structure (as long as it has tetragonal 
symmetry) and the prefactor /i only depends on the an- 
gular slope of the gap function and Fermi velocity at the 
node. Numerical evaluation shows that the free energy 
( p^ ) gives rise to a centered rectangular FL structure, 
aligned with x or y axes, but now with the apex angle 
P < 60°, depending non-monotonically on the field. This 
suggests that there will be an additional transition at low 
temperature related to the non-analytic behavior of the 
response tensor. 

In order to map out the complete equilibrium H-T 
phase diagram we have carried out a numerical compu- 
tation of the FL structure using the full expression for the 
response tensor Q(k), as given by Eqs. (||) and ([T^). We 
find that free energy has two local minima for centered 
rectangular lattices aligned with two high symmetry di- 
rections shown in Fig. [l|, as expected from the tetragonal 
symmetry of the problem. Which of the two becomes 
the global minimum depends on temperature and field. 
The results are summarized in Fig. ^ For high temper- 
atures the exact result agrees well with the one obtained 
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from the long wavelength free energy (||). The defor- 
mation of the lattice from perfect triangular grows with 
decreasing temperature, reaches a maximum, and then 
falls. Maximum distortion occurs around t ~ 0.3, attain- 
ing /3 ~ 70° at lOT. Extrapolating this field dependence 
(see inset to Fig. ||) , the FL should become square around 
H « SOT, but this field is outside the domain of validity 
of the London model. At lower temperatures the distor- 
tion decreases but instead of going all the way back to 
triangular at t* , the lattice undergoes a first order phase 
transition to another centered rectangular lattice rotated 
by 45° and with /3 < 60°. Further decrease of tempera- 
ture causes the angle to grow again. We note that precise 
temperature at which it crosses 60° depends on field, but 
for all fields is close to t* — 0.19, as predicted by the long 
wavelength approximation. At yet lower temperature we 
predict another first order transition to a centered rect- 
angular lattice along x (or y) with /3 < 60°. The free 
energy difference between the two minima is very small 
in the region where the 45° rotated lattice wins. It is thus 
likely that real system will remain in the metastable state 
and the experiment would detect only a smooth crossover 
from a lattice with /? > 60° to the one with /3 < 60°. 

The present calculation can be easily generalized to 
treat the effects of Fermi surface anisotropy. As men- 
tioned above tetragonal anisotropy will not modify the 
T universal behavior but may lead to quanti- 
tative changes at higher temperatures. Orthorhombic 
anisotropy, on the other hand, will modify even the 
T limit. We expect that it will, to leading order, 
merely rescale the coordinate axes, leading to the same 
structures as described above stretched by the appropri- 
ate factor It may further remove the degeneracy 
between two equivalent lattices related by 90° rotation. 
Another source of anisotropy neglected in our calculation 
is the non-linear Meissner effect studied by Yip and Sauls 
associated with the shift of quasiparticle spectrum 
due to the supcrflow. Within the quasiclassical picture 
this would contribute terms ~ {\dx' B\'^ + \dy' B\'^) to the 
London free energy at T = 0, where x' and y' are 45° 
rotated coordinates. Our initial numerical results ||l9|] in- 
dicate that while not completely negligible, these terms 
will not significantly alter the behavior of FL reported in 
Fig. ||. Also neglected in our calculation is the effect of 
electronic disorder, which will remove the non-analyticity 
of Q(k) at longest wavelengths, just as small finite tem- 
perature would. Since the FL is most sensitive to Q(k) 
at finite k ^ (/is the vortex spacing), we expect our 
predictions to be robust with respect to weak disorder. 

In conclusion we have described distinctive features of 
the vortex lattice geometry associated with unconven- 
tional pairing and non-local response. Near our pre- 
dictions are consistent with the existing phenomenologi- 
cal work while at low T we predict novel effects 

intimately related to the nodal structure of the order pa- 
rameter, which are not contained in GL-type theories. 
These will occur simultaneously with other unique effects 



predicted previously, such as the ^ y H dependence of 
specific heat [|o). Existing experiments probing the FL 
geometry in cuprates |^,^ provide a somewhat conflicting 
picture and their theoretical analysis is complicated by 
the orthorhombic anisotropy and presence of twin bound- 
aries. We hope that the present theory will encourage 
more systematic experimental work, preferably on un- 
twinned or tetragonal materials. 

After this work was completed we learned about a 
preprint by Kosztin and Legett which discusses ef- 
fects of non-locality on the effective penetration depth 
in d-wave superconductors. Where overlap exists their 
results appear consistent with ours. 

The authors are indebted to A. J. Berlinsky and Z. 
Tesanovic for helpful discussions. This research was sup- 
ported by NSERC, the CIAR and NSF grant DMR- 
9415549 (M.F). 



[9: 

[lo; 

[11 
[12; 

[is: 

[14 

[is: 
[16: 

[17 

[is: 

[19 
[20 
[21 



A. A. Abrikosov, Sov. Phys. JETP 5, 1174 (1957). 

For review see: Anisotropy Effects in Superconductors, 
ed. by H. Weber (Plenum, New York, 1977) 
K. Takanaka, Ref. §. 

G. Blatter et al, Rev. Mod. Phys. 66, 1125 (1994). 

B. Keimer, et al, Phys. Rev. Lett. 73, 3459 (1994). 

I. Maggio-Aprile, et al., Phys. Rev. Lett. 75, 2754 (1995). 

K. Takanaka and K. Kuboya, Phys. Rev. Lett. 75, 323 

(1995); M. Ichioka, et al, Phys. Rev. B 53, 2233 (1996). 

A. J. Berhnsky, et al. Phys. Rev. Lett. 75, 2200 (1995); 

M. Franz, et al, Phys. Rev. B 53, 5795 (1996); J.-H. Xu, 

Y. Ren and C.-S. Ting, ihid, R2991 (1996). 

E. H. Brandt, Phys. Rev. Lett. 78, 2208 (1997). 

I. Affleck, M. Franz and M. H. S. Amin, Phys. Rev. B 

55, R704 (1996). 

M. R. Eskildsen, et al, Phys. Rev. Lett. 78, 1968 (1997). 
V. G. Kogan et al, Phys. Rev. B 55, R8693 (1997). 
See e.g. M. Tinkham, Introduction to Superconductivity 
(Krieger, Malabar, 1975). 

E. M. Lifshits and L. P. Pitaevskii, Statistical physics II 

(Pergamon Press, New York 1980) pp. 208-216. 

V. G. Kogan et al, Phys. Rev. B 54, 12386 (1996). 

E. H. Brandt, J. Low. Temp. Phys. 24, 709 (1977); 73, 

355,(1988); and Phys. Rev. B 37, 2349 (1988). 

D. J. Scalapino, Phys. Rep. 250, 329 (1995). 

S. K. Yip and J. A. Sauls, Phys. Rev. Lett. 69, 2264 

(1992). 

M. H. S. Amin, I. Affleck and M. Franz (in preparation). 
G. E. Volovik, Sov. Phys. JET P 58, 469 (1993). 
I. Kosztin and A. J. Leggett, (|;ond-mat/9702199|). 



4 



